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Abstract 



We calculated some of the critical exponents of the directed percolation univer- 
sality class through exact numerical diagonalisations of the master operator 
of the one-dimensional basic contact process. Perusal of the power method 
together with finite-size scaling allowed us to achieve a high degree of ac- 
curacy in our estimates with relatively little computational effort. A simple 
reasoning leading to the appropriate choice of the microscopic time scale for 
time-dependent simulations of Markov chains within the so called quantum 
chain formulation is discussed. Our approach is applicable to any stochastic 

process with a finite number of absorbing states. 
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In its original formulation |]T], |2|, the so-called directed percolation (DP) conjecture stated 
that all continuous phase transitions about a single absorbing state in single-component 
systems with a scalar order parameter are in the DP universality class of critical behaviour 
@. In this form the conjecture has been confirmed in a host of model systems, amongst others 
the basic contact process j|, [5], Schlogl's models for autocatalytic chemical reactions ]I|, 0, §|, 
and a phenomenological classical field theory of high energy hadronic collision processes 
|7|, |, H [10| . Further investigation revealed that the DP universality class is robust enough 
to accommodate more general models, with more than a single component 
as well as with multiple, in some cases infinitely many absorbing states 



ni pi, O, P| 



15, 16, 17 . Even 



some nonequilibrium growth models without absorbing states were found to share some of 
the DP exponents 



21, 22 



[lq , |19| , P0|| . For recent reviews see 
The basic contact process (CP) [|], |] may be viewed as a model for the spread of an 
epidemic amongst individuals living in a (^-dimensional orchard. In this process, a healthy 
individual becomes infected at a rate proportional to the number of its infected neighbours, 
whilst an infected individual X becomes healthy at unit rate. Pictorially, in one dimension 
it is defined by the elementary processes X00 ^ XX0, 00X -4 0XX, X0X Q XXX, 
and X — > 0. As A increases from zero, the basic CP suffers an extinction-survival phase 
transition in all dimensions, the upper critical dimension being d* = 4. There is not an 
exact evaluation of the critical points A* for d < d* to date, but there are some narrow 
bounds: in one dimension it is known that 1.539 < A* 



23 



< 1.942 

In this work we were concerned with the accurate determination of the critical point 
and some of the critical exponents of the one-dimensional basic CP through exact numerical 
diagonalisations of its master operator. Our method is based on the standard matrix power 
method |24j applied to a discrete-time version of the continuous-time Markov chain, taking 



advantage of the presence of an absorbing state. Combined with finite-size scaling |23| and 
modern extrapolation techniques [26|, the method allowed for a high degree of numerical ac- 



curacy within quite reasonable computational efforts. Successful attempts at the application 
of phenomenological renormalisation group ideas to directed problems on the lattice dates 
back at least to the work of Kinzel and Yeomans || , and have from time to time reappeared 
in the literature with varied levels of sophistication [^7]. The application of the 'quantum 



chain formulation' of Markov chains on the lattice |$| to study time-dependent properties, 

In this work we provide a simple 



29 



however, remained scarce; for a recent example, see 

reasoning leading to the appropriate choice of the microscopic time scale for simulations that 

should be of value to anyone interested in similar calculations. 

The starting point of our approach is the master equation on the lattice. Let AcZbe 
a one-dimensional lattice of |A| = L sites with periodic boundary conditions. To each site 
£ G A we attach a random variable ne taking values in a finite set uj = {0, 1, . . . , N — 1} cN, 
the state of the whole lattice being given by n = (m, 712, . . . , til) G Q = to . In the basic CP 
sites can only be healthy or infected, thence N = 2. Given positive real numbers r(n, n) 
denoting the rates at which the collision n — > n occurs, the master equation governing the 
time evolution of the probability P(n, t) of realisation of the configuration n at instant t 



reads 

^P(n, t) = J2 ( r ( n > fi ) P ( fi > *) " T(n, n)P(n, *)]. (1) 

n 

We now introduce vector spaces in the description of the above equation. To do this we turn 
u into C^ and write 

|P(t)) = ^P(n,t)|n) (2) 

n 

for the generating vector of the probabilities P(n, t) = (n|P(t)), with {|n)} the orthonormal 
basis diagonal in the occupation number representation. We are in this way providing the 
space of generating functions with a Hilbert space structure. We then rewrite Eq. ([J) as 

A|P(t)) = -H\P(t)) (3) 

where H is given by 

tf = ££>(n,n)|n><n| (4) 

n n 

with H(h, n) = — r(n, n) and H(n, n) = X^n^n r(n, n). H is but the infinitesimal generator 
of the Markov semigroup U(t) = exp(-tH) of the continuous time Markov chain {n(t),t > 0} 
defined by the set of rates r(fi, n). The spectrum of H lies in the complex right half-plane, 
and since if is a real matrix, its eigenvalues are either real or come in complex conjugate 
pairs. The steady state of H has eigenvalue zero, and for finite systems it is, up to symmetry 
degeneracies, unique. If we realise the algebra of operators of Q in terms of spin S Pauli 



matrices [28], the master operator H of the basic CP can be seen to be equivalent to a spin 
S — |, three-body non-Hermitian quantum chain. 

The lowest gap E^ 1 ' — E^°> = E^ in the spectrum of H may be used to perform a finite- 



size scaling analysis in the same way as one is used to do in equilibrium problems ||25|| . Let 
us briefly review some formulae. Around the critical point A > A*, the correlation lengths of 
the infinite system behave like 

£H oc ^ oc (A - A*)- 1 -" oc (A - X*)- U±z (5) 

where £ii and £j_ are the correlation lengths respectively in the time and space directions, 
i/ii and uj_ are the corresponding critical exponents, and z = u»/u± is the dynamic critical 
exponent. Notice that in the interacting particle systems literature it is more usual to find z 
defined as z = 2u±/vn. For finite systems of size L, according to the usual finite-size scaling 
assumptions [^] we expect 

ql = L-^(\X-Xl\L 1 ^ (6) 

where zl and u±,l are the finite versions of z and u± and $(x) is a scaling function with 
$(x ^> 1) ~ x u u. On general grounds one expects lim^oo A^, Zl, v±,l = X*,z, v±. From 
Eqs. (H) and (0) we obtain 



ln(L/P) ln(L'VL) ZL [ ' 



which through the comparison of three different system sizes L' < L < L" furnishes simul- 
taneously X* L and zl. According to Eq. (|), the exponent z/j_ jL may be calculated through 
v ~^l = z l + \{Cv,l + Cl,l"), where 



In 

Cm.n = — 



d & 9X ) J { d ^,M/ dX 



ln(N/M) 

with the derivatives calculated at A = X* L . Of course £,7^ = Re{E ( L }. 

Several possibilities exist to proceed with the calculation of the gaps. When if is a 
symmetric matrix, Lanczos diagonalisation becomes the method of choice. For stochastic 
processes, however, H is in general non-symmetric, and although there exist non-symmetric 
variations of the Lanczos algorithm, they are either much more memory demanding or are 
intrinsically unstable [23], Since we are interested in only one very precise eigenvalue, we 



choose to work with the power method, which requires only matrix-by-vector multiplications 
that can be carried out using high precision data types. 

In order to apply the power method, we first define the matrix T = 1 — tH . This 
matrix should be viewed as the time evolution operator of a discrete-time Markov chain 
whose transitions take place at intervals r. For T to be a stochastic matrix, its elements 
ought to satisfy < T(n, n) < 1 and ^~T(n, n) = 1. This last condition is always 
satisfied, simply because Y2n H(ii, n) = 0. The first condition, however, demands that r _1 > 
max{|if (ri, n)|} = maxjif (n, n)}. Since we will calculate the eigenvalues of T by an iterative 
procedure, convergence will occur at a maximal rate if we choose r as large as possible, and 
we take r _1 = 1.01 x maxjif (n, n)}. It is not advisable to take r _1 = max{if(n, n)} 
because this will zero some elements in the diagonal of T, thus introducing cycles in an 
otherwise acyclic Markov chain. The spectrum of T lies in the unit circle, with the steady 
state corresponding to the eigenvalue one. Notice that the spectra of H and T appear in 
reverse order. The above choice of the microscopic time scale r for the discrete-time Markov 
chain is equivalent to the requirement that the probability of two transitions taking place in 
time t is negligible, of o(t). This makes the approximation U(r) = exp(—rH) ~ 1 — tH to 
the Markov semigroup exact in what regards the dynamics, i.e., it preserves stochasticity, 
and the time-evolved vectors \P(t + r)) = exp(-rH)\P(t)) and \P(t + r)) = (l-rH)\P(t)) 
coincide up to o{r). 

The construction of the matrix T presented above together with the appropriate choices 
of r with the purpose of studying time-dependent properties of stochastic processes is well 
known in queueing theory, where with some minor refinements it is called the method of 
uniformisation |30] . 



In the basic CP in finite volume, the steady state is given by the unique absorbing state 
|0) = |0,0,...,0). This a priori knowledge of the steady state of the process makes it 
possible to calculate the second largest eigenvalue of T, that corresponding to the gap of H, 
by simply orthogonalising the iterated vector of the power method against the steady state 
at every iteration. We thus end up with an implementation of the power method that reads 

|P(m+l)) = (T-|0)(0|)|P(m)> (9) 



with m in units of r. We expect that after enough successive applications of the above 
relation it reaches a fixed point E* \ P*) =T\P*) where E* = (P*\T\P*)/(P*\P*) = l-r£ (1) 
and \P*) is a linear combination of one-particle states. The advantage of being dealing with 
an absorbing state in contrast to, e.g., a numerically determined steady state, is that the 
former is usually a 'pure' state, as is our case, or a rather simple combination of states with 
known coefficients, e.g., white noise, a fact that both minimises the inevitable numerical 
round-off errors and saves computer time, turning the calculations more reliable and fast. 
In our calculations, we consider an eigenvalue to have converged if it coincides more than 64 
times with its predecessors in more than one part in 2 112 ~ 5.2 x 10 33 . 

The complete characterisation of the DP universality class requires, besides z and u±, one 
more exponent to be calculated, the other exponents following from well known hyperscaling 
relations |9|, [L^, ^T[. One calculable exponent is 8, defined through the asymptotic behaviour 
of the survival probability at the critical point A* as 

PsurvW = Ys P ( n > *) = l ~ P ( ' f ) * ^ ^^ + ^^ + •••)• ( 10 ) 

A logarithmic plot of P S mv{t) versus t for an infinite system should be a straight line at large 
t, the slope of which is 8. The correction exponent 8' seems to be given by 8' = 1 [[H], [31]. 
For finite systems, however, the spectrum always has a gap, and the survival probability 
ultimately enters a regime of exponential decay ruled by the finite gap. In order to follow 
the time evolution of the process before it gets trapped into an absorbing state, we take a 
small t in the definition of T and successively calculate \P(m + 1)) = T\P(m)). We choose 
t _1 = 1000 x max{iJ(n, n)}. The initial state was given by |P(0)) = |1, 0, . . . , 0), the state 
with a single particle located at the origin. A plot of P SU rv(t) versus t for a system of L = 16 
sites is shown in Fig. |l| From that figure we clearly see that after an initial transient in which 
the 'highest modes' are washed out, the curve enters a regime of almost pure algebraic decay 
until the gap in the spectrum manifests itself and we begin to observe the expected late times 
exponential behaviour. This is more clearly appreciated in the inset of Fig. [I], which shows 
the derivative of logP surv (t) with respect to logt, that is the instantaneous value of Sl- 

Our results for A*, z, u± and 8 are summarised in Table [TJ. The L = oo values in 
this table were obtained through the Bulirsch-Stoer (BST) extrapolation scheme p6| , with 



ujbst the free parameter of the algorithm chosen so as to minimise the difference between 
the penultimate entries in the BST extrapolation tableaux. The derivatives in Eq. ^ were 
calculated with a nine points symmetric difference formula with an 0(h 9 ) error using steps of 
h = 10~ 9 [0, whilst the values of 8 were obtained through linear regression fits to log -P S urv(/) 
in the region of algebraic decay of P SU rv(^), typically through 100 data points with points 
separated by a time interval At = lOOr. The numbers in Table [I| confirm that the basic CP 
belongs to the DP universality class of critical behaviour with a high degree of accuracy. The 
uncertainties associated with the extrapolated numbers are mainly due to finite-size effects 
and corrections to scaling, as well as to the extrapolation procedure itself. Our numbers 
compare well with those obtained by other means, namely, high-temperature expansions on 
a closely related 'reggeon quantum spin chain' model 0, other finite-size scaling studies 
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Figure 1: Survival probability at A = A^ for a system of L 
instantaneous value of the critical exponent Sl- 
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16 sites. The inset shows the 



time-dependent operator perturbation calculations P3[ , and extensive series expansions and 
differential approximants analyses [Q. The best values of z, u±, and 5 to date are given by 
z = 1.580 745(10), u ± = 1.096 854(4), and 5 = 0.159 464(6) El. Our estimate of A* is as 



precise as those available in the current literature on the basic CP [23, 33]. Although our 
determination of 5 is not as precise as that of z and u±, it nevertheless is precise enough 
to discriminate amongst the universality classes that are likely to arise in systems with 
absorbing states pl[. An alternative to estimate 5 is to determine each 8l as the value of 
the derivative of the corresponding logP surv (t) with respect to logt at its inflexion point, 
see Fig. [I]. This provides a more 'local' determination of 5l than that obtained by fitting a 
straight line to logP surv (t) x logt over tens or hundreds of points, and may improve the final 
result once the uncertainties are correctly assessed. We intend to pursue this alternative in 
a future, more thorough study. 

A remark about symmetries. In the calculation of the gaps of H, one can take advan- 
tage of any symmetries of the process to reduce H to block-diagonal form. The numerical 



Table 1: Finite-size data and extrapolated values for the the critical point A* and the expo- 
nents z = un/u±, v±_ and 8 of the one- dimensional basic CP. The numbers between paren- 
theses represent the estimated errors in the last digit of the data, whilst those data without 
an associated error are numerically precise to the figures shown. 



L', L, L" 


K 


z L 


V±,L 


h 


7,8,9 


1.629 092 086131 


1.495 084128194 


0.963 208 351697 


0.1657(2) 


8,9,10 


1.632 522 345 029 


1.502 980 235 818 


0.977844 866 308 


0.1656(2) 


9,10,11 


1.635178 201359 


1.509 743 775 238 


0.989427140 315 


0.1654(2) 


10,11,12 


1.637262 542 035 


1.515 558 577190 


0.998 833 401056 


0.1653(2) 


11,12,13 


1.638 921714266 


1.520 588 023 063 


1.006 632 974979 


0.1651(2) 


12,13,14 


1.640 260 494445 


1.524 967656 357 


1.013211430 501 


0.1649(2) 


13,14,15 


1.641354409414 


1.528 807324384 


1.018 839 364 776 


0.1647(1) 


14,15,16 


1.642 258 557889 


1.532195 497121 


1.023 712 402 567 


0.1645(1) 


15,16,17 


1.643 013 687274 


1.535 203 516105 


1.027975 558429 


0.1644(1) 


16,17,18 


1.643 650 350 303 


1.537889180610 


1.031738 664206 


0.1643(1) 


17,18,19 


1.644191762 995 


1.540 299 611765 


1.035 086 466 706 


0.1641(1) 


18,19,20 


1.644655 789 991 


1.542 473 490 560 


1.038 085 430 711 


0.1640(1) 


L = oo 


1.648 96(2) 


1.580 77(2) 


1.09681(2) 


0.162(2) 


[^bst] 


[1.071] 


[1.171] ' 


[0.895] 


[2.701] 



diagonalisation is then performed in the sector of lowest gap with a considerable economy 
of computation. Internal symmetries, like U(l) and Z(N) symmetries, usually separate the 
dynamics into sectors corresponding to the closed classes of the stochastic process, with 
each block a stochastic transition matrix governing the dynamics within the given sector. 
Geometric symmetries, however, like translation and reflexion symmetries, generally lead to 
block operators that are not stochastic due to the occurrence of 'artificial' combinations of 
states. If one decides to make use of geometric symmetries, care should then be exercised in 
properly weighting the basis vectors in order to interpret the resulting symmetry-invariant 
\P{t)) as a vector of probabilities. In this work we explored the translational invariance of 
the basic CP on a periodic lattice in order to achieve a reduction of order 1/L on the size of 
the matrices we needed to diagonalise. 

In summary, we have successfully applied the power method to the one-dimensional basic 
CP, and obtained accurate values for the critical point A* and the critical exponents z and u±, 
together with a good estimate of the critical exponent 5. The method is fast, yields accurate 
estimates for the critical point and some of the exponents, and is easily coded. Given that it 
took less than 500 h of CPU time (running at 300 MHz) to complete Table [I], and that in the 
largest cases it took less than 21 Mb of memory to conduct the calculations, the method seems 
to be very competitive. Extension to processes with more than one absorbing state as well 
as in more than one dimension is immediate. It would be of interest to refine the calculations 
of 5 as well as to try to calculate other exponents by the same methods. In particular, it 
would be very interesting to establish relationships between dynamical exponents like 8 and 



the spectrum of H . We are at the moment pursuing these objectives, and intend to release 
our results soon. 
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